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CHAPTER 1 : INTRODUCTION 

Scope 

Potential flow methods have been used successfully for 
the last two decades in the preliminary design of partial or 
complete aircraft configurations. Predictions by numerical 
schemes based on potential flow analyses include such aero- 
dynamic characteristics as wing load distributions, surface 
pressure distributions, engine duct flows, and some stability 
derivatives , among others . 

Results of such computations concerning surfaces affect- 
ed by trailing vortex sheets of lifting surfaces are correct 
only for simple configurations where the effects of the loca- 
tion and shape of the vortex sheet are secondary. The shape 
of the vortex sheet is usually assumed to be flat. Realisti- 
cally, however, its configuration changes continuously in the 
downstream direction at least until roll-up of the sheet into 
concentrated vortex cores is complete. The process of roll-up 
is rather complex. Mutual interaction among the elements of 
the sheet depends on their relative positions, however, the 
configuration of the sheet is unknown prior to the complete 
solution including its effect. 

Much of the change in the sheet shape occurs within a 
downstream distance from the generating wing equal to half 
the wingspan. The purpose of this study is to develop a 
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numerical technique of modeling the vortex sheet with a 
deformable surface definition, along which a continuous vor- 
tex strength distribution in the spanwise direction is ap- 
plied, so that by repeatedly modifying its shape, its true 
configuration is approached, in the proximity of its generat- 
ing wing. 

Design problems requiring the inclusion of a realistic 
configuration of the vortex sheet are numerous. Some 
examples are discussed in the following. 

Control effectiveness and stability derivatives 

In the early stages of aircraft design, horizontal and 
vertical stabilizers must be sized fairly accurately to en- 
sure aircraft controllability within the flight envelope. 
Downwash and sidewash angles at zero angles of attack and 
sideslip and their rates of change with these angles are pre- 
dicted in practice using empirical relations (reference 11 
and 31, for example) which are not based on general aircraft 
configurations, but rather on crude parameters such as wing 
sweep and dihedral angles, aspect ratio, etc. The real 
governing factor is the wing loading distribution and changes 
in it with Mach Number and attitude with respect to free 
stream direction. This distribution results in a free vor- 
tex sheet extending in theory from the trailing edge of the 
generating wing to an infinite distance downstream of the 
wing. The vorticity is constant streamwise, and varies with 
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chordwise distance. Roll-up of the sheet occurs about areas 
of concentrated vorticity called the vortex cores (4, 23) 
and if in the vicinity of the empennage, such cores will 
affect surface pressure distributions. Two cases are 
represented. 

Longitudinal stability As illustrated in Figure 1, 
the deflected flaps of the aircraft result in a redistribu- 
tion of lift on the wing, and thus a concentrated vortex 
core pair emanates at the outboard edge of the flaps that 
drastically changes the downwash at the horizontal stabilizer. 
Shock induced separation of the outer wing panel in transonic 
aircraft causes a similar effect, which aggravates the insta- 
bility due to the forward displacement of the center of lift 
of a sweptback wing. 

Lateral stability Figure 2 depicts a yawed twin- 
engine aircraft. The dip in the wing lift distribution due 
to the nacelles, caused by inviscid as well as viscous 
effects, leads to a pair of concentrated vortex cores which 
are not symmetric due to yaw. The sidewash angle distribu- 
tion along the stabilizer is altered. In addition, the left 
side of the rear fuselage is closer to center of the inboard 
left vortex than the right side is to the center of the in- 
board right vortex. A pressure differential and thus a 
destabilizing yawing moment result. 
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Figure 2. Lateral destabilization due to nacelle side 
forces and complex wing loading 
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Canards 

An effect similar to the above should be investigated 
for canard-configured aircraft. The vorticity from the 
canard could affect the lift distribution of the main lift- 
ing wing, and subsequently alter the aerodynamic performance 
of the wing under certain conditions. 

Propellers and helicopters rotors 

Good aerodynamic design for propeller or rotor after- 
bodies (nacelles and other solid boundary surfaces) requires 
simultaneous consideration of all components. The effect of 
these lifting surfaces is felt on other solid boundaries 
through the vortex sheets they shed. Therefore, a more de- 
tailed knowledge of the shape of the vortex sheet and hence 
the induced flow field is necessary. For these cases, a 
quasi-static analysis is required. Such studies are also 
necessary for designs which must result in low aerodynamic 
noise. 

Trailing vortex hazard 

Although this area has been studied extensively recently 
(7, 9, 23), the present method could be extended to its treat- 
ment, insofar as determining the locations and intensities 
of each of the trailing vortex cores. The interest in this 
case is primarily in the far downstream region and to promote 
early dissipation of vortex energy. Farther downstream, 
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viscous effects become more pronounced, and since the 
current method is inviscid, vortex merging and dissipation 
must be treated by some other means. Nevertheless, it is 
necessary to establish the initial roll-up which is essen- 
tially inviscid, to enable the treatment of merging and 
dissipation. 


Background 

The present method is based upon earlier work in three 
major areas of computational methods, all of which have 
undergone considerable development in the last two decades. 

A brief historical review is presented in this section. 

Panel methods 

In the early sixties, a computational method was devel- 
oped to predict nonlifting potential flow about arbitrary 
three-dimensional bodies by placing a source distribution on 
the outer surface of the body. The source strength is approx- 
imated by piecewise constant strength over flat "panels" 
which approximate the body surface (15). At about the same 
time, vortex lattice methods were being developed (26) for 
very thin lifting surfaces. The approach is similar except 
that the camber surface of a lifting wing is represented and 
a discrete horseshoe vortex lattice is used. Later, the two 
techniques were combined to solve arbitrary three-dimensional 
lifting potential flow configuration problems. Examples of 
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work published in this area include those of Rubbert and 
Saaris (28). The solution of the flow is obtained by solv- 
ing for the singularity (source and doublet) strength dis- 
tributions on the entire surface. The problem is reduced 
to a set of simultaneous linear equations which can be solved 
non- trivially by imposing Neumann boundary conditions at the 
solid boundaries and the Kutta condition at the lifting sur- 
face trailing edges. The induced effects on the surface or 
elsewhere can then be calculated. More recent developments 
were directed to higher order singularity (source and doublet) 
distributions, in order to improve the predictions and also to 
facilitate surface paneling (20) . Excessive care for the 
panel arrangement is thus no longer required. 

In the present study, the doublet (or vortex) strength 
distributions are known (obtained by some panel technique) 
for the test wing. The spanwise distribution remains the 
same along the sheet at any streamwise station, downstream 
of the trailing edge. The same type of doublet distribution 
panels are used to determine induced effects to obtain the 
rolled-up shape, except that these panels are constantly 
reshaped and refitted. 

Vortex sheet roll-up 

A comprehensive survey of computational methods for 
lift-generated wakes is given by Rossow (23). The vortex 
sheet (or wake) trailing behind the generating wing is divid- 
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ed in three regions: roll-up, plateau, and decay. In this 

study, the roll-up region is the one of interest. The prob- 
lem of roll-up was addressed by Betz (4), then by Westwater 
(32), where approximations were introduced to eliminate the 
unnecessary complications of viscous effects. Many techniques 
have been developed since to overcome the shortcomings using 
artificial means. The major problem was the use of discrete 
vortices to approximate the vortex sheet, rather than a con- 
tinuous distribution as the present method proposes. 

The interest in obtaining the rolled-up shape of the 
wake was motivated by the advent of large aircraft. The 
wake left by a Boeing 747 for example, during climb or 
descent is extremely hazardous for following aircraft and 
extends for several kilometers. Many of the methods develo- 
ped have been intended for finding means of alleviating that 
hazard, such as by modifying the wing loading of the generat- 
ing aircraft to prevent or reduce vortex merging (part of the 
roll-up mechanism for complex wing loading) and thus speed up 
the decay process. 

Parametric bicubic surface representation 

Representation of general curved surfaces was developed 
for computer aided design of aircraft (13) and automobiles 
(5). It is necessary in this study to use curved panels, and 
it was found that the Ferguson patches are best suited for 
modeling the vortex sheet surface. 
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General Description of the Proposed Technique 

Since the purpose of this method is to establish the 
shape of the trailing vortex sheet in the vicinity of the 
wing, it is assumed that the strength of the vortex sheet 
is known. Curved and deformable panels of known piecewise 
linear vorticity strength distribution (equivalent to quad- 
ratic doublet distribution) are patched together to simulate 
the vortex sheet, (Figure 3). The geometry of these panels 

is modeled by using parametric bicubic patches with second 

2 

derivative continuity (C ) across their boundaries. 

Initially, the surface is flat. The induced velocity at 
each node in the network is calculated, and some nodes can 
be displaced accordingly, resulting in a new shape, with the 
remaining nodes following. The procedure is repeated with 
the new shape until the entire surface has been relaxed to 
the final shape. (A complete description is in Chapter A.) 

These induced velocities can include the effects of the 
bound and free doublet panels as well as any solid boundary 
whose consideration is desirable, such as a fuselage, trail- 
ing lifting surface, nacelle, etc. However, the solution for 
the entire flow must be re-computed at each relaxation step. 
In this study, only the doublet panels of the trailing and 
bound vorticity are considered, since the inclusion of other 
bodies is merely a programming and computing effort. In 
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Node 1,1 



Node 1,2 


Boundary 

v=l 


u and v are the parameters 
in the bicubic expression 
for physical coordinates. 


Figure 4. Panel nomenclature 
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addition, the roll-up effect on the lift distribution on the 
generating wing itself is assumed to be negligible (30) . 

Flow Idealizations Summary 

The following assumptions are made in this study: 

Potential flow The flow is everywhere irrota- 
tional and incompressible except at the boundary surfaces 
(specifically the vortex sheet). 

Symmetry of the generating wing and vortex sheet 
Sideslip could easily be included for yawed configurations. 

Inviscid flow The viscous effects become pro- 
nounced downstream of the area of concern. 

Tip flow No tip flow for the wing is considered, 
since the wing character itself is not of primary concern, . 
although its loading near the tip will be somewhat affected. 

Steady state Oscillatory and other time-depen- 
dent flows must be treated quasi-statically if this technique 
is to be used. 

Wing lift distribution The wing lift distribu- 
tion is based on the assumption of a flat trailing vortex 
sheet. As indicated above, the roll-up effect is negligible 
at least for the purpose of this study. In Rossow (23), 
noticeable effects are shown which could be included with 
considerably more computational effort. 
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CHAPTER II: VORTEX SHEET GEOMETRY 


The vortex sheet is mathematically considered to be a 
surface of finite width extending from the lifting surface 
to infinity, across which there is a jimp in the potential 
function <J> . This jump in <p is equal to the doublet strength. 
The surface must be represented by a mathematical model which 
provides for at least slope continuity and allows for the 
curled or spiral shape of the rolled-up sheet including infi- 
nite slope tangents, as well as slope constraints at the 
edges of the sheet. 


Parametric Bicubic Representation 


The position vector x=(x,y,z) of a point P on the 
surface within one patch is given by the interpolant 

x = I I a, u K v (1) 

k=0 £=0 KIC 


where u and v are the independent parameters , such that 
0<u, v< 1 , and a^ are the polynomial coefficients. Note 
that the vector coefficient a^ consists of three elements, 
one for each of the x, y, and z coordinates. The sixteen 
vector coefficients can be determined by specifying x,x u> x v 
and x uv at four comer points or nodes. (The subscripts u,v 
refer to the partial derivatives and , respectively.) 


9v 
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The entire surface is defined by a rectangular array of 
patches. Every four adjacent patches share one common inter- 
ior comer point or interior node. Around the periphery, 
every two adjacent patches share one common boundary point or 
boundary node , except that the comer nodes belong only to 
the comer patches. This array of patches will be referred 
to as a network . An illustration is shown in Figure 5. 

In Figure 5(a), the geometry for a network of m x n 
patches is shown in the parametric space. A total of 
(m + 1) (n + 1) nodes result. The interior nodes are those 

with index i = 2,3,...,m and j = 2,3 n. The others are 

boundary nodes. Each parameter changes from 0 to 1 within 
a patch, and the coordinates of the parametric space shown, 
u' and v 1 , are the cumulative values of the u, v. The bound- 
aries are lines of constant u or v parameters. Normal to the 
u = constant boundaries, are the tangent vectors x^ for each 
node; and normal to the v = constant boundaries, are the 
tangent vectors x v , parallel to the v = constant and 
u = constant lines respectively. These tangent vectors 
required must be specified according to considerations to be 
discussed in Chapter 4. 

The transformation in equation (1) takes a point in the 
parametric space to a corresponding point in the physical 
space. It is necessary to determine the remaining 2m + 2n -4 
tangent vectors at the boundary nodes, plus 2(m-l)(n-l) 
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tangent vectors at the interior nodes. This tangent vector 
information, together with node coordinates, is sufficient 
to enable the computation of the sixteen vector coefficients, 


[a, ] (in matrix form) according to the matrix relations 


U kl ] -QXQ' 


where X = 

x 00 X 01 x v00 x v01 

x 10 X 11 x vl0 x vll 

x u00 x u01 x uv00 x uv01 

_ x ul0 x ull x uvlO x uvll 

and Q = 1 0 0 0 

0 0 1 0 

-3 3 -2 -1 

2-2 1 1 

The matrices for the y and z coordinates are similar (12) . 

The letter subscripts indicate partial differentiation, 

and the numerical ones indicate the value of the parameter 

at the patch nodes: the first for u and the second for v. 

3 x 

For instance, x u q^ refers to (u=0, v=l) . 


Advantages of parametric bicubic patch panels 

The main advantage of a higher order patch geometry for 
an aerodynamic singularity panel resides in its curvature. 

A much larger number of flat panels is required to obtain the 
same error in position. For instance, the cubic representa- 
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tion of a 90° circular arc results in a radial error of +0.13% 
of the true circle using one curve segment (12). To obtain 
the same accuracy with straight segments, eleven would be 
required. 

The cubic requirement is necessary to allow for second 
derivative continuity and inflection points (12) which are 
inherent characteristics of the nature of the vortex sheet. 

In addition, the parametric bicubic is widely used in air- 
craft design, (3, 18), and its technology is well-established. 

The fact that the vortex sheet surface changes directions 
on itself in the spanwise direction to result in the spiral- 
like shape, renders explicit functions cumbersome to use. 
Implicit functions are equally cumbersome since they would 
require solutions of nonlinear equations at each point. The 
parametric representation eliminates these difficulties, and 
moreover handles vertical slopes without special provisions. 

Computations of the interior tangent vectors 

With reference to Figure 5(b), consider the set of nodes 

(i,j) with a fixed value of i, i* , that is, the nodes (i',1), 

(i',2), ..., (i',n+l). The patch edges passing through these 

nodes constitute a space curve, which is a composite para- 

2 

metric cubic which must possess the C property, if the entire 
2 

surface is C . In the case of any i' curve, we obtain a curve 
whose parameter is v. Similarly, another set of curves, j', 
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are obtained by varying the parameter u. By differentiating 
the equation for a curve with constant v: 


k 

x(u) = E a, u 

( 3 ) 

k=0 k 


x (u) =, ^-,ka,u k ^ 
u k=l k 

( 4 ) 

3 

(u) = Z k(k-l)a,u k ^ 

1 o K 

( 5 ) 


The requirement is that at the junction of the ith and (i+l)th 
segments : 


x i (l) 

- x i+ \ 

0) 


^(l) 

u 

= x i+1 
u 

(0) 

(6) 

x 1 (1) 
uu 

->i+l 

= X 

uu 

(0) 



Equation (6) must be satisfied to achieve second order continu- 
Now, as for the patches in Equation ( 2 ), the curve 
polynomial coefficients can be written in terms of the 
position and tangent vectors at the nodes: 
gent vectors at the nodes: 


a 0 r 0 
a l ■ r u0 

a 2 = 3(r l- r 0> ' 2r u0 - r ul 
a 3 = ' 2<r 0 - r l> + r u 0 + r u 


( 7 ) 


Substituting these values in equation ( 5 ), then in the third 
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of Equation (6) 


6x 1 - 6 x 1 + 2X 1 + 4X 1 = 6x i+1 - 6x i+1 - 4$ i+1 - 2x i+1 
t>x Q dx 1 + zx uQ + 4x ul ox 1 ox Q 4 x uQ zx u1 

But, due to the first two of Equation (6), 


x 1 + 4x i+1 + $ i+1 = 
uO + ^ x u0 + x ul J * x l O' 


Now, Xq refers to node i, Xq = x i refers to node i + 1 and 
# | *1 

refers to node i + 2. Thus for nodes 2, 3, ..., m, the 
system of equations results: 


J 1 + 4S i+1 + J i+2 = ScJ 1 * 2 -?-) 
u u u 


x 1,2,..., in™ 1 


( 8 ) 


where i here denotes the ith node rather than patch as done 
earlier. 

If x* and x™ + ^ are known, then the remaining m-1 tangent 
vectors at the interior nodes can be found efficiently by 
solving the above tridiagonal system by the Thomas' algorithm, 
( 1 ). 

The above is repeated for all n+1 curves, then, in the 
transverse direction, for all m+1 curves in a similar fashion. 

The remaining derivatives x uv , also called the cross 
derivatives, or twist vectors, have zero values along lines 
of minimum or maximum surface curvature. For simplicity, 
this assumption was made and, as indicated by (12) , results 
in a negligible amount of error. 
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In the same reference, it is shown that the patches are 
2 

continuous in C across their common boundaries if the above 
curve continuity is achieved. 

Geometric Properties of Interest 

The main purpose of a patch is to serve as a surface 
position interpolant, i.e., given a parameter pair of values, 
a corresponding set of physical coordinates of a point on the 
surface of patch i is obtained: 

x = ( x j y , z ) = I Z af.u^v 5, 
k=0 £=0 

such that 0 £ u, v £ 1. 

In addition, a tangent vector to the surface at the 

point is a linear combination of a pair of tangent vectors , 

_ -*■ 

x = -r— and x, T = ■*— which are obtained merely by differentiat- 
ing with respect to the parameters u or v. If x y ^ ax.^, 
where a is a scalar, then the vector product x u * x v is a 
vector normal to the surface at the given point. 

Elemental area As illustrated in Figure 6, the 
elemental area dA is given by: 

dA = | x u du| • | x v dv | sin 0 = |x u * x v |dudv 

= l i <y u z v ■ y v z u ) + j(x v z u ' x u z v } + k(x u y v ■ x v y u )l dudv 

A A A 

where i,j,k are the unit vectors in the x,y,z directions. 
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Integration of a function over the surface of a patch 

It will be necessary to integrate the induced effects of 
the vorticity which is distributed over the patch area, and 
expressed as a function f of u and v. To evaluate the 
integral 

I -// f(u,v)dA, (10) 

S 

where S is the boundary surface bound in one patch, and dA is 
the elemental area, the relation above is used, thus 

U 

I =//f(u,v)|x x x Idudv. 

00 u v 


(ID 
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CHAPTER III: SINGULARITY SOLUTIONS FOR POTENTIAL FLOW 

Governing Equations and Solutions 


For steady, irrotational motion of an inviscid perfect 
fluid, the velocity perturbation potential function satisfies 
the Prandtl-Glauert equation: 

(1-M 2 ) d> +4 + <j> =0 

v <» / Y xx y yy y zz 

which is valid for small perturbations, but not for transonic 
flow. For subsonic flow, by replacing x by x(l-M ) 2 , the 
equation reduces to Laplace’s equation: 


v^<j> - $ + 

r T xx 


4 > + 

yy 


zz 


0 


This equation also holds for incompressible irrotational flow 
in general. 

Applying Green’s theorem, it can be shown (2) that 
Laplace's equation can be converted to the integral equation: 

♦p-fc"t-M4 + *-|n‘5” ds (12 

where: P is a point in the flow field domain V, 

S is the boundary of V, 

n is the unit vector normal to S at a point Q on the 
surface, and 

r is the distance from P to Q. (Refer to Figure 7.) 
This equation is a singularity solution to Laplace's equation. 
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The two terms in the integrand contain 

Twr and 'k ( ^7 ) 

which are the potential functions' of a point source and a 
doublet, respectively. The terms <j> and represent the 
source and doublet strengths at point Q, thus equation (12) 
can be rewritten as: 


* U> Sf + y(X) 3n ( ||.-| >]dS 

where: £ is the position vector of field point P, 


x is the position vector of boundary point Q, 


o(x) source strength distribution on the boundary S, 
and y(x) doublet strength distribution on the boundary S. 

In panel methods, a and y are approximated by piecewise 
constant, linear or quadratic distributions over each surface 


finite element, called a panel, so that for K panels: 


♦<*> "jHr + •*! Is7 ldsl (U) 

where is the area of the ith panel. By properly specify- 
ing boundary conditions at a number of points (K) , called 
control points, selected in this case on the surface, a set 
of simultaneous linear equations can be obtained and solved 
for o i ,y i - 

Usually, source panels are used for closed surfaces and 
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solid boundaries, while doublet panels are used for thin 
lifting surfaces, or for lift generation on solid boundaries 
as well as for barriers such as vortex sheets to render V 
simply connected. 

In this study, only the doublet solution will be of 
concern, even though source distributions may be included to 
represent solid boundaries and their effect on the roll-up 
process of the vortex sheet. 

Higher Order Doublet Panels 

Unlike the panel method techniques where the objective 
is to solve for the singularity strengths, in the present 
case, the wing lift distribution, a function of the spanwise 
location, is known a pA-ioA.1. It could be determined for a 
given configuration from some panel method computer code, 
experimentally, or by some other means. The problem here is 
to determine the induced effects of the vorticity at various 
key points of the flow field, in order to determine displace- 
ments and thus modify the shape of the vortex sheet. Constant 
doublet strength panels are equivalent to discrete ring vor- 
tices, and a substantial approximation results. In the 
interest of maintaining adequate accuracy without an excessive 
number of panels of constant doublet strength, a bilinear or 
biquadratic doublet distribution is used. 
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Bilinear: 

y T (u,v) = y, + y 9 u + vi oV + y,uv 
1 l Li 4 (15) 

Biquadratic : 

2 2 2 2 2 2 
( u »v) = UjCu.v) + y 5 u + y g v + UyU v + Uguv + yqU v 

» 

where u and v are the arc length in the strearawise and span- 
wise directions, respectively, of the vortex sheet. 

In the free portion of the vortex sheet, the doublet 
strength does not vary strearawise, and u is a constant there. 

Equivalence of doublet and vortex panels 

The induced velocity at a point P due to an elemental 
vortex panel of area dS is obtained from the Biot-Savart law: 

= -T^rYdS (16) 

v I ir 

where: r is the relative position vector of P with respect 

to the panel centroid. 

A 

t is a unit vector in the direction of the vorticity. 

Y is the local vorticity strength. 

The velocity potential at a point P due to a doublet 

A 

panel of local strength y with normal vector n is: 




A A 

It can be shown that y*t is equivalent to n x Vy , (16), 


thus producing the same potential function in V. Both 
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approaches are used. In reference (27) , the doublet panel 
is preferred because the doublet strength is a scalar. 

However, a biquadratic approximation is used to reduce the 
approximation error, since the first order error term is 
nonexistent, resulting in a more complex panel model. 

Reference (18) suggests the use of a bilinear vorticity 
distribution, which is an equivalent approximation to the 
biquadratic doublet panel. 

To assess the approximation error, consider a two- 
dimensional flat vortex sheet as in Figure 8. The induced 
vertical velocity at point x on the sheet segment from x^=0 
to x^=c is given by: 

, c y (x, )dx-, 

V z< x - 0) * K ' --- - - X - < 17 > 

Expanding y(x^) in Taylors' series about x: 

y(x^) = y (x) + •g--(x) • (x^-x) + ^-^(x) (x-^-x) 2 + 

Substituting in Equation . (17) : 

V x>0) = h [ + ^ (x) + h ^ (x)(x r x > + o 2 <x 1 -x)]dx 1 

"k^ (x) '^ + ai< x >- c+ i$< x > l < x i-x>“*i> + ^° 2 < x f x >- 
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If point x is the midpoint of a panel, the first and third 

terms in the square brackets vanish. Thus, in order not to 

dv 

exceed second order error, the derivative ^ must be nonzero. 
This is achieved by using a bilinear y distribution. Biquadra- 
tic distribution will not reduce the error. In addition, since 
the bilinear distribution over each panel provides (continu- 
ity of value) for y, there will be no concentrated vorticity 
at the panel edges. 

Analysis similar to the above can be used to demonstrate 
that source panels possess the same behavior, and bilinear 
source distributions would result in a second order error in 
the distance from the control point. 

Derivation of the vortex panel induced velocity for the free 
vorticity 

Considering the geometry of one panel represented by a 
parametric bicubic, the vorticity vectors are tangential to 
lines of constant values of the parameter v, which are space 
curves lying in the vortex sheet, oriented in the general 
streamwise directions as shown in Figure 9. By virtue of 
Helmholtz' first vortex theorem, the strength will not vary 
along these curves. Therefore, each v-line is a vortex 
filament of infinitesimal strength y(y-^)dy^, where y-^ is the 
y-coordinate of the generating wing trailing edge. As indi- 
cated earlier, y = y(y^) is known. For simplicity, it will 

* 

be approximated by a piecewise linear distribution in v 



Figure 9. Vortex panel configuration for vorticity 
along a v=v, curve 
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between the nodes of the free vortex sheet's leading edge, 
which coincides with the wing's trailing edge. For the 
panels with index j , 


Y (v) = y(y\j_^) + [y(yj) - (18) 

at any point on such panel whose position vector x is given by 
Eq. (1). The unit tangent vector is given by 

x u 

t = — (19) 



The vorticity induced velocity at a point whose position vector 
is f due to the vorticity in panel i,j is obtained by integrat- 
ing Eq. (16) using Eq. (11): 


= 


11 




v t r 7 


( 20 ) 


with m x n is the number of patches representing the sheet in 
the stream- and spanwise directions, respectively. 


Derivation of the vortex panel induced velocity for the bound 
vorticity 

The formulation is identical to that of the free vorticity 

/\ 

with the exception that the vector t does not coincide with 
x y , but is a linear combination of x^ and x v » Equation (20) 
is thus applicable to the component along inducing V^. 

The other component along x v induces a velocity vector V 
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11 

ff 

00 


x u x (£-x) 

l* u l |t-x | 3 


Y^u.v) |x u x 


x | dudv 


2 


1 

fn 


11 x x (|-x) 

"o |S V | ||-x| 3 ^2^ u,v ^ x u x x v |dudv 


also , 


and for the entire sheet consisting of mxn panels: 


m n . . 

^ = Z Z $ 1,J 

i-1 j =1 


(21a) 


(21b) 


( 22 ) 
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CHAPTER IV: THE VORTEX SHEET 

Review of Methods of Rollup Prediction 


Betz 1 Theory 

The earliest attempt to analyze the behavior of vortex 
systems in a manner applicable to the prediction of the then 
known phenomenon of vortex sheet rollup was introduced by 
Betz in 1932 (A). The theory is based on the Kutta-Joukowsky 
theorem of vortex lift on a body L, about which a circulation 
T exists: 

L = pvr 

where p is the fluid density, and V is the linear velocity of 
the body relative to the fluid perpendicular to L. Thus, if 
a system of such bodies interact, with the bodies becoming 
infinitely small, and if there are no solid outer boundaries 
to exert a force, the net lift forces and moments of these 
interacting vortices must be zero. As a consequence, and as 
a consequence of Helmholz' theorem, the following conservation 
rules can be written to relate the vorticity in the sheet as 
it leaves the wing trailing edge, with that of the rolled-up 
sheet: 

1. Vorticity: 


b/2 

/ 

y 


dr (y) 

dy 


dy 


R 

/ 


dT(r) 


dr 


(23) 
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The left-hand side of the equation refers to the 
initial flat sheet at the wing trailing edge, y 
being any spanwise station and b the wing span. 
The right-hand side refers to the rolled-up vor- 
tex core, r being a radial position in the vortex 
core from its center and R the core radius. 

2. The first moment of vorticity: 


b/2 

-/ 

y 


dr (y) 
dy 


R 

ydy = f 
r 


dr (r) 
dr 


rdr 


(24) 


As a consequence, the centroid of the vorticity 
located at : 


y 


1 

r, 


b/2 

/ 


dr (y) 

~3y 


ydy 


where 


F 0 = 


b/2 

/ 

y 


dr (y) 
dy 


dy 


remains at a constant spanwise station. The vertical 
position, however, will change since there is an 
external force, the wing lift, applied to the fluid 
elements which results in a downward motion of the 
centroid. 

3. The second moment of vorticity about the centroid: 



40 


b / 2 tI 21 <y-y) 2 d y - / r 2 dr (25) 

y y r 

These relations can be used to determine the structure of the 
rolled-up vortex core, that is, the position of the core 
radius r in which a portion of the sheet from a station y to 
the tip is contained. Rossow (25) shows that 

r = y - y 

Wings with nonsimple loading, such as for flaps extended 
configurations can also be treated using a method developed 
in reference (8) . 

Unfortunately , these rules and methods are helpful only 
for fully developed vortex cores, since the time history of 
the rollup cannot be predicted. 

Time History of Rollup 

In 1935, Westwater (32) developed a technique for obtain- 
ing the vortex sheet shape by representing it with discrete 
vortices of finite strength and infinite length, and computing 
their positions at successive time increments. These positions 
are equivalent to flow streamlines and the time increments 
correspond to streamwise stations. Computations are performed 
in the Trefftz Plane, that is, the infinite filaments are 
assumed straight and perpendicular to the plane. The bound 
vorticity is ignored and thus the problem is reduced to two- 
dimensions.. The induced velocity at a point where any one 
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vortex filament intersects the Trefftz Plane is the sum of 
the induced velocity vectors of each vortex filament except 
the one in question. The point is moved by an incremental 
displacement equal to the time increment times the induced 
velocity vector. This is done for each filament and the 
process repeated for various streamwise stations to obtain 
the time history of the configuration. Plots of Westwater's 
results are given in his paper as well as in (30) and (9) 
for an elliptic loading using 20 filaments for the entire 
span. 

Discretization seems to lead to several problems. The 
elimination of the filament at which the induced effects are 
being computed is both valid and necessary for a finite or 
discrete system, however, it is a poor representation of a 
continuous system since it creates a "gap" in the sheet. 

Should the number of filaments be increased, the two vortices 
adjacent to the point in question would become too close and 
result in excessive induced velocities. Another problem is 
the increased strength of the tip filaments, compared to the 
remaining ones. Moore has found, as reported by El-Ramly (9), 
that such tip filaments will circle around each other, imply- 
ing that the sheet intersects itself. 

Some of these idealizations were removed by Butter and 
Hancock (6). Three-dimensional effects were considered. A 
line vortex at the quarter chord was introduced to represent 
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bound vorticity, and the trailing vorticity is only semi- 
infinite. A similar effort was conducted by Hackett and 
Evans (14). The discretization effects are the same as 
above. 

Reference (9) notes the work of Neilsen and Schwind, 
where some treatment is done to alleviate the discretization 
effects, namely, when two vortices come too close to each 
other, they are combined into one at their centroid. 

Other researchers whose works are surveyed by El-Ramly 
(9), concentrate on fully rolled up vortex cores, the dis- 
tances for full roll-up and other related parameters not 
directly related to the present study, for both inviscid and 
viscous solutions. The inclusion of viscous effects to pre- 
dict decay of vorticity and the merging of co-rotating cores 
has received considerable attention in the '70s decade when 
the wake hazard of large aircraft became manifest. Again, 
these works are surveyed in (9, 23). Experimental work was 
conducted at Iowa State University (reference 19 ) , to study 
the merging of vortex cores which are fully developed. There, 
two separate wings were used to generate a pair of rolled-up 
vortex cores , in order to simulate multiple cores generated 
by a single wing of a large aircraft during climb or approach 
near airports, where the wake hazard is serious. The merging 
of such cores may be desirable since it can lead to vortex 
cores with diffused vorticity representing a reduced hazard 
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to trailing aircraft. Much of the research (9, 24) deals 
with the injection of auxiliary vortices which merge with 
the primary tip vortex, as well as methods of assuring 
speedy decay. For a specific aircraft configuration, it 
is helpful to know where such cores will be situated upon 
completion of rollup, and the relative magnitudes of their 
strengths. Then, analytical or experimental predictions of 
merging and decay can be simplified. 

The present method in essence is similar to that of 
Butter and Hancock (6) , but removes many restricitons which 
lead to over-simplifications. Mainly, a continuous distribu- 
tion of vorticity, truly three-dimensional effects of the 
vortex sheet shape, and a lifting surface are used with the 
help of the building blocks described in Chapters II and III. 
Furthermore, the technique is compatible with the advanced 
aerodynamic paneling techniques in use today, allowing ready 
incorporation of rollup effects, if desired. 

The remainder of this chapter describes the method used 
from the physical and computational standpoint. 

Sheet Vorticity Strength Distribution 
Panel Vorticity 

At any point P on the vortex sheet (see Figure 10), the 
vorticity strength vector y can be resolved in two components 
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8 r 8 r 

Y-, = -5 — and y 0 = -r— — , which are tangential to the 
i- dS^ Z dS2 

u=constant and v=constant curves at point P, respectively. 

T is the circulation in the sheet at point P, and s^ and S2 
are arc lengths in the direction of and respectively. 
Note that y^ and are not necessarily perpendicular to each 
other. It is necessary to express these vectors in terms of 
their values at the patch nodes, using a bilinear distribu- 
tion as indicated earlier. For the magnitude, 

Y (u, v) - U 21 Y i+l,j + U 12 Y i,j+l + U 22 Y i+l,j+l (26) 

where, 

y . . are the magnitudes of y-, or Yo at nodes (i,j), 

U^ = l-u-v + uv, 

U 2 = u - uv , 

ti^2 = u - uv, and 

U 22 = uv * 

A A 

For the direction, the unit tangent vectors t-^, t2 are 
merely the tangents to the v=constant and u=constant curves, 
respectively, so that: 

t l = X u / I x u I > and 

t 2 = V> x vl 


( 27 ) 
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Computation of nodal vorticity strengths 

It is assumed here that: (1) the wing spanwise loading 

and (2) the chordwise loading distributions at various span 

locations are known. A number of wing stations (span loca- 

dr 

tions) are selected, and the values of ^ at these sta- 

tions obtained, where V is the bound vorticity strength at 
the spanwise station y. The value of T is obtained from 
elementary vortex lift theory (Kutta-Joukowsky) : 


r - 7 V ~ C C* 


where C and are the local chord length and lift coeffi- 
cients, respectively. 

Between a pair of consecutive wing stations, j and j+1, 
y-^ is approximated by a linear distribution in y: 


Y 

1 


(y) = 


<^i ) i+i - <2lh 


y. 


j+i 


- 7: 


(y - yj) + (Yx)j 


The change in bound vorticity strength in this interval, 

T, +1 - Tj must be equal to the strength of the shed or trail- 
ing vorticity. As a direct result of Helmholtz' first 
vorticity theorem: 


y j+i 

r j+i - r j ‘ L Y i dy 
y ± 


y±+i - y 


5 h t'n’j+l + 1 
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Tj is distributed in the chordwise direction in a fashion 
similar to that of the load distribution, , where Cp 

is the pressure coefficient, and the subscripts refer to lower 
or upper surfaces of the wing. Figure 11 depicts typical dis- 
tributions for two categories of airfoils. In addition, the 
following relation holds: 

= / Yodx 
J 0 z 


at any wing station j . y ^ denotes the component of vorticity 
proceeding in the spanwise direction. It is again approxi- 
mated by a piecewise linear distribution (Figure 12), i.e., 
linear over a chord segment from point i to point i+1. 


Y 

2 


(x) 


(y 2 } i-H 
x i+l ' x 


(Y 2 >i 


(x - x.) + (Y 2 >i 


Between two consecutive stations j , j+1, a linear blending 
function is used to interpolate, and thus the bilinear dis- 
tribution is used. The same interpolation method is used 
for Y-, . 

C 


Now, 


since T 


/ y^dx, the piecewise interpolant 
0 Z 


results in a trapezoidal rule summation: 


r . 

j 


™ x i+i - x i , , , 

2 C (Y 2 5 i+1 + (r 2 , l I 
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Figure 11a. Typical chordwise vortex strength 
distributions 



Figure lib. Approximation of chordwise loading 
distribution using a linear 
interpolant 
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Figure 12. Typical spanwise load distribution for a wing 
and piecewise linear approximation of corre- 
sponding shed vorticity 
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Since the distribution is known, the (Y 2 ^ can t> e expressed 
as multiples of jr,, the maximum Y 2 va lue along the chord, 
then 


r j = 7 y 2 . + ( Vi ]4x i = f ' Y 2 

■J 1 = 1 


(28) 


or, 



(29) 


where : 



l=l 




and 


Ax. = 


X 


i+1 


x . 


- V t 2 


(30) 


Therefore, the vorticity flow through the panel edge (1), on 
Figure 13, has the circulation 


r l = ‘< Y 2>i+l + <)'2 > i 1 j 4x i 


“ Y 2 [(Y 2>i+l + (Y 2 ) i 1 j Ax i 


(31) 


and through panel edge (2) 


r 2 = [(v 2 ) i+l + (Y 2 ) i 1 j+l 4x i = y 2 [(y 2 ) 1+1 + <Y 2 ) i , j 4x i <32) 
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Similarly, through panel edge (3), the vorticity flow has the 
circulation : 


r 3 " t'Vj+l + Ay j 

and through panel edge (4): 


(33) 


r 4 = t(y i>j+i + 

The following conditions must be noted: 


l. iy + r 3 = r 2 + iy 

(35) 

2- Y 1 = Y 2 = 0 

(36) 


along the leading edge, and 

3. Yl - 0 (37) 

along the centerline or plane of symmetry. 

If T is known at the centerline, and y-^ is known at all 
wing stations, which are determined a. pnioni. along with the 
chordwise loading, then by virtue of Equations 28-30, 36, 

7*2 at each station, and thus Yi at all nodes can be obtained. 

By virtue of Equations 31-36, the values of y^ at all nodes 
can be obtained too. 

In the free portion of the vortex sheet, the wake y ^ ~ 
everywhere, and the values of y-^ remain unchanged from those 
at the trailing edge. 

In the parametric patch representation scheme, all 
variables are expressed in terms of the parameters u and v. 


0 



53 


which are used to represent the streamwise and spanwise 
directions with proper transformations, respectively. 

As a matter of general interest, in the panel methods 
where the circulation is everywhere unknown, it is customary 
to express the nodal circulation values in terms of a subset 
of values, usually one per spanwise station. The induced 
effects are then added to the source panels representing the 
solid boundari'es . An additional number of equations equal 
to the number of unknowns (in this case the yJ, values) is 
required. These are obtained by specifying boundary con- 
ditions at an equal number of points, for example, where the 
Kutta condition may be enforced. The problem of the vector 
nature of vorticity is thus eliminated. 

Remark on spanwise positioning of nodes 

Parametric bicubic patches of the Ferguson type used in 
this study guarantee second derivative continuity. However, 
the node spacing has an important effect on the quality of 
fit. If two neighboring patches have widely varying lengths, 
the tangent vectors at the common boundary will be too high 
for the smaller one, causing loops, and resulting in excessive 
surface areas as well as area related integrals over the 
patch. 

On the other hand, the piecewise linear distribution of 
vorticity dictates another criterion in the selection of nodal 
positions. A linear segment should be terminated before the 
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error in the linear representation is excessive, in order to 
maintain the higher order terms at a minimum. 

Induced Velocities 

It is necessary now to compute the induced velocity due 
to the entire vortex sheet, bound and free portions, at a 
number of points on the free portion, namely, the patch nodes. 
Equation 22 is used. The induced velocities of each panel, 
Equation 21 (a generalization of Equation 20) , are summed for 
all panels representing the vortex sheet. If other singulari- 
ties are present, such as solid boundaries represented by 
source panels or other lifting surfaces represented by vortex 
panels, their induced velocities would be added, too. Source 
panels induced velocities are given by an expression similar 
to Equation 20, except that the vortex strength is replaced 
by source strength, and the tx term is dropped (18). 

The surface integral in those equations is not suitable 
for closed form evaluation, unfortunately, and a numerical 
quadrature is employed (Chapter V) . 

Symmetry 

Computations can be greatly reduced by taking advantage 
of the symmetry of the configuration about the y=0 plane. 
Letting the physical quantities be expressed for the "left- 
hand half" of the wing in terms of those on the "right-hand 
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half" (primes are used for the left side) , then: 

x ' = x 

y' = -y 

z ' = z 

for the vorticity unit vectors: 

^ /V a- 

t, = x, i + y j + z k 
1 u J u J u 

A AAA 

t,' = -x i + y j - z k 

1 u J u J u 

A AAA 

t 0 =xi+yj +zk 

2 v J v J v 

A A A. A 

ti = -x i + y i j - z k 

2 v J v J v 


(39) 


( 40 ) 


for the distance from the elemental vortex to a point in space: 


6x ' = 6x = E, - x 

6y' =n+y; Sy=n-y (41) 

6z ' = 6z = ^ - z 

Again, % is the position vector of the point at which the 
induced velocity is to be computed, x is that of the vortex 
element, and subscripts u and v denote partial derivatives 
with respect to these parameters. In Figure 14, these various 
quantities are depicted. 

The displacements in the wake elements are mirror imaged 
on the left side from those computed on the right side; thus, 
the panel symmetry relations (38) hold. 



-a 


Plane of 
Symmetry 



Figure 14. Symmetry of geometric parameters 
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Procedure to Compute Rollup 

Initial conditions and sheet length 

The proposed method is an iterative procedure where an 
initial shape for the free portion of the vortex sheet must 
be assumed. The most logical one would be to extend the trail 
ing edge in the direction of the free stream velocity, since 
all sections of the sheet has gone through the trailing edge 
as they were being shed from the wing. The sheet is extended 
to a downstream . station far enough so that the induced effects 
at the areas of interest (e.g. , the empennage) become negli- 
gible. An additional extension is required to promote adequat 
influence for the rollup of original portion. In the latter 
portion, accurate rollup is not of essence, rather, its Simula 
tion of the fact that the vortex sheet is in reality semi- 
infinite in length is the intended purpose (Figure 15) . 

The free portion of the sheet is divided in three regions 

- Region I in which accurate estimates of induced velocities 
are needed. 

- Region II is an extension to such distance beyond which the 
rollup effects are negligible everywhere in region I. 

- Region III is a further extension to such distance beyond 
which vorticity effects are negligible in region I and II. 

The various streamwise stations (or spanwise or u=constant 
curves) are those along which the patch corners or nodes are 
located, and are numbered from 1 at the trailing edge to 




Figure 15. Breakdown of trailing vortex sheet for computation purposes 
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Computation scheme 

The following steps are performed to obtain an approxima- 
tion of the rolled-up shape of the vortex sheet. 

1. The induced effects of the entire initial sheet at 
the nodes of stations 1 and 2 are computed. 

2. The induced velocities at station 1 nodes are re- 
quired for computing the downwash and sidewash angles 
at the nodes, such that the resultant velocity 


V- (V co + AV X U + AV V j + AV Z k 

3 X j 7 j Z j 


(42) 


where the V' s are the induced velocity components , 
coincides with the tangent vector x . The magnitude 

, j 

of x^ should be approximately equal to the arc length 
of a v=constant curve through the particular node on 
station 1 between this node and the corresponding 
one on station 2. Thus: 




(43) 


where £ is approximated by the chord length l x 2j -x ljl 
for the node pair in question. The nodal displace- 
ments are zero. 

3. At the nodes of station 2, displacements at the nodes 
are computed. The time elapsed from shedding to the 
position along station 2 is: 
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try . = (x 9 . - x. . ) /V j = 1,2,..., n+1 (44) 

During this time t, the induced velocities vary 
from those induced at station 1 to those at station 
2. The average velocity is used, so that the dis- 
placement for each node at station 2 is: 


: 2,J< V i.J +V 2,j ) 


j = 1,2,. 


(45) 


Thus, the position vectors are obtained for station 


2 nodes: 


x 2,j “ x 2,j + Ax 2,j 


The superscript refers to the computation cycle, zero 
being the initial condition. The treatment of the 
tangent vectors at the three remaining boundaries is 
presented later in this chapter. 

4. For all the remaining stations, 3 to m^^, the incre- 
mental nodal position vectors are set equal to those 

I 

at station 2. The reasoning is that the vortex sheet 
past station 2 has passed through that curve at some 
earlier time. 

5. The second computation cycle is similar to the first, 
except that now the third station (i=3) is modified 
according to the following: 
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A^ is an incremental velocity at station 3 to repre- 
sent the acceleration of the particles between sta- 
tions 2 and 3. Now, the time elapsed is 


c i,j = (x i,j - x i-i,j )/v = 


(48) 


Then, the incremental displacement Ax is: 




t. .(AV 2 . . + AV 2 . , .) 

i.J i.J l-l. J 


(49) 


6. ' For stations 4 to the incremental nodal posi- 

tion vectors are equal to those of station 3 nodes, 
and so are the tangent vectors at the tip. 

7. The cycle is repeated m.^ times, (for i values up 
to and including m^) , to obtain the relaxed wake 
shape, so that for the kth cycle: 

A^k . = ^ _ ^k-1 

i,J i.J i.J 


H.j - (x i,j - x i-i,j )/v . 


->-k -»-k+l . 1 . , 

x — x + t t • • ( AV 


+ 4 V i-i,: ) (50> 


i = 2,3,4,..., k+1 


Particle motion due to a vortex is not linear. An 
assumption of linearity is valid only for small time segments. 
In order to relieve this restriction, circular motion about 
each elemental vorticity is approximated as follows. 

In Figure 17, let point P be a field point under the 
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influence of a vortex element at point Q. The distance 
between them is p. The true updated position of P is P’, 
resulting from circular motion about Q. P^ is the new 
position using linear motion approximation which can be 
easily obtained, but is obviously inadequate unless p > > | 6 V • 1 1 . 
Point ?2 is obtained by multiplying the vector - § by the 
ratio p/|?^ - §| . Then, P" is obtained by multiplying the 
vector ?2 ” ^ by the ratio |?2 * ?| » so that ?" - ? 

has the same magnitude as the induced effect 6^*t. For 
moderate angles (up to tt/ 2), P" is a reasonable approximation 
of P'. 

Tangent vectors 

The free portion of the vortex sheet is refitted after 
each displacement. An evaluation of the tangent vectors in 
transverse direction to the boundary is required. In step 2, 
the tangents x y were evaluated at the wing trailing edge. 

The tangents x^ at station m^^ (the last or pseudo infinity 
station) are assumed parallel to the free stream velocity. 

Due to symmetry, the tangents x v at the symmetry plane are 
perpendicular to that plane. The tangent vectors of the 
remaining boundary, the tip vortex, are discussed next. 

Tangent vectors at the tip vortex 

First, the direction of the vector is determined as 
follows: By referring to 'Figure 18a, consider the current 




Figure 18a. Exact tangent vector computation 




/ 


/ Approximate 
tangent vector 
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section in the sheet represented by curve C. The tip point A 
has the position vector x = (y, z). The y and z coordinates 
are those of interest since changes in the x-direction are 
assumed negligible. Point B is a point on curve C just in- 
board of A whose position vector is (x - dx) = (y-dy, z-dz). 
The parameter v assumes the value v and v - dv at points A 
and B, respectively. Curve C' represents the updated shape 
of the vortex sheet after the next cycle, that is, after a 
time t has elapsed since the particles have moved from the 
upstream station, as in Equation 48. The induced velocity 
vectors at points A and B are V and V - d^, leading to new 
positions A' and B' , with position vectors x' and x 1 - dx' . 
Thus, the tangent to the curve at A' has the components dy' 
and dz ' , where 

dy' = dy + dV^t and 

dz ' = dz + dV t , 
z 

or the components y^. and z^, the length components of the 
tangent vector, where 

X = v(y v + (V y ) v t and 

z v " v(z v + (V z>v t) * 

v is a scale factor which can be obtained by equating the 

2 2 k 

magnitude [y^. + j to arc length of the last segment 

of the curve, approximated here by its chord length 

t( yn+l - >n )2 + <z n+l ‘ Z A )2]% - Theref ° re . 
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.2. <y;+l - yn )2 + <Z A+1 - 2 n> 2 

v ■ 2 ! 7 

(y + (V ) t + (z + (V ) t) 

7 v y v v z V 

Note that 

3V 3V 

<V v - a ? 2 and <Vv - -at 

can be obtained by differentiating Equation 22 with respect 
to v using Leibniz' rule. Since the limits of integration 
are constant, such differentiation is merely performed on 
the integrand. 

The evaluation of would be performed by numerical 
quadrature in the same fashion as ^ itself, necessitating 
the near doubling of computations to be performed. Thus, an 
approximation is adopted as shown in Figure 18b. A point (B) 
is selected inboard of the tip (A) at a parameter value Sv 
less than that at the tip. Sv is chosen small so the chord 
B'A' represents the tangent vector of the updated curve C' 
at the tip A' . The position of B' is computed in the same 
manner as the panel nodes. 
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CHAPTER V: 

DESCRIPTION OF COMPUTER PROGRAM 
Introduction 

In this chapter, the computer program to perform the 
computations required for the algorithm described in Chapter 4 
is presented. It is written in FORTRAN 77 in a structured, 
modular form for easy expansion, upgrading, and adaptation to 
an aerodynamic analysis system. 

Major Module Flowchart 

A flowchart showing the major functions is shown in 
Figure 19. Each major function or module is detailed further. 
Three sets of streamwise panel rows are used: 

- The bound vortex panels which are not updated. 

- The updatable panels at which the displacements are 
calculated. 

- The slaved panels which are updated to the last station 
evaluated. (Correspond to those of Region III. See 
Figure 16.) 

The "updatable row" repetition constitutes the relaxa- 
tion cycle iterations. The entire wake update is repeated 
for the rows of Regions I and II which the induced velocities, 
displacements and tangent vectors are re-evaluated. The 
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Figure 19. Major function flowchart 
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results could be output on an external file in the same format 
as the inputs and the program re-run using those for a better 
approximation, although this was not done in the present study. 
As will be shown in Chapter VI, the results obtained from the 
first iteration were quite satisfactory. In the following, 
the various details are discussed. 

Surface fit 

This module merely fits parametric cubic curves in the 
u and v directions independently, resulting in the tangent 
vectors at the nodes in both these directions. The twist 
vectors are not obtained, for simplicity, and due to their 
minimal effect on the quality of fit. However, if necessary, 
they should be computed here. The Thomas' algorithm is used 
to solve the set of tridiagonal equations as explained in 
Chapter II. 

Quadrature points 

This module interpolates the wake surface within a 
patch for a set of 25 points, (see Table below), used for 
quadrature, to obtain the position and tangent vectors, 
x, x u and x v , Then, . the latter two are used to obtain the 
normal vectors x u * x v and their magnitudes, and store these 
data for use by other portions of the program within the 
patch loop. This module uses another one to obtain the poly- 
nomial coefficients of the patch, given the nodal positions 
and tangents. The latter are merely the partial derivatives 
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of the position polynomials, with respect to the parameters 
u and v. These data are used in obtaining the deformed or 
stretched patch area in order to compensate for the stretch- 
ing, as well as for evaluating the induced velocity integrals. 

Numerical integration 

The Gaussian quadrature is used. The integrand is 
evaluated at twenty-five points in the interior of each panel 
forming a grid of five points in u-direction along five con- 
stant v curves. The values of u and v are shown in the table 
below. The elemental area is approximated by a weighting 
factor for the integrands a^a^ , the values of a^ corresponding 
to the grid point is also shown in the Table. The integral 

11 

// f(u,v) dudv 
00 

is approximated by the double summation 

5 5 

£' E f(u. , v.) a. a. 
j=l i=l 1 J 1 J 


Gaussian Quadrature Parameters 


u. or v. 



1 

2 

3 

4 

5 


0.046910077 

0.230765345 

0.5 

0.769234655 

0.953089923 


0.11846344 

0.23931434 

0.28444444 

0.23931434 

0.11846344 
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Induced velocities 

For each of the quadrature points , the values of the 
vorticity are interpolated, and the relative position vector 
(to the node at which the velocity is being calculated) evalu- 
ated for one half of the symmetrical vortex sheet, as well as 
the other opposite half. The induced velocity vector induced 
by the current patch at the entire set of nodes is obtained 
by use of the Gaussian quadrature. These are stored and incre- 
mented for all patches of the network. This is done so that 
the patch quadrature points are evaluated only once per relaxa- 
tion cycle, to reduce the amount of computations. In this 
fashion, the bulk of the computations is done in the induced 
velocity quadrature. 

Typical run statistics 

One test case used 15 x 7 panels. The 15 streamwise 
rows comprise 3 bound, 7 slaved, and 5 updatable sets. The 
run required 105 CPU seconds on the Iowa State University 
NAS-AS/6 computer, constituting most of the thirteen dollar 
charge. Should this method be used in conjunction with a 
potential flow paneling program, two iterations are to be 
sufficient, add a small percentage to the total cost. 
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CHAPTER VI: 

RESULTS AND DISCUSSION 

The results of the computational method of the present 
study are presented in two parts: 

1. A wing planform tested for vortex sheet shape 
visualization is modeled numerically and compared 
with the test results to demonstrate the validity 
of the method. 

2. The span loading of the wing is altered arbitrarily 
to simulate deployed flaps (or stalled outboard 
panel) is modeled and the results presented to 
demonstrate the ability of the system to handle 
complex loadings. 

The first part uses the results of a wind tunnel test con- 
ducted in the Boeing Co. Research Wind Tunnel, which is 
described next. 


The Test Wing 


Description of the test 

The planform of the wing is shown in Figure 20, along 
with the spanwise load distribution, measured experimentally 
using pressure taps at various stations. The load distribu- 
tion used in the numerical model must be a piecewise linear 
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distribution of bound vorticity. This modeling is described 
in the following section. 

The flow visualization technique used for this experi- 
ment consisted of injecting steam and liquid nitrogen in the 
airstream at a location upstream of the wing tip. The nitro- 
gen flow rate is adjusted to bring the air and water vapor 
mixture close to the dew point, so that a small drop in tem- 
perature will cause condensation of the vapor and result in a 
mechanical mixture of air and water droplets. This effect 
will take place if an adiabatic pressure drop occurs, specif- 
ically along the vortex sheet, where the induced velocities 
are locally high. If the wind tunnel is of the open type 
and the atmospheric relative humidity is high, and the air 
temperature is very low, liquid nitrogen would be unnecessary. 

Light is applied to a section of the stream past the wing 
trailing edge through a narrow slit so that only the water 
droplets at the particular section are illuminated and can 
be photographed. 

A complete description of this and other methods used 
for this program can be found in references (21, 29). The 
only available photograph from the test was taken at a 
section 1.25 spans behind the wing's trailing edge, and is 
shown in (29), along with computational results discussed 
later. The flow visualization is illustrated in this report 
as the dashed line in Figure 22f. 



Twist angle, degrees 
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Figure 20a. Test wing planform and twist 



Figure 20b. Test wing span loading (experimental) 


203 . 2 mm 
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Numerical modeling of wing loading 

The necessary inputs to the computer program are piece- 
wise linear vorticities for each panel of the vortex sheet. 
In this section, the bound (or fixed) portion thereof is 
discussed. 

Since no chordwise loading was published in (29), a 
distribution similar to aft loaded airfoils was assumed. 
Three chordwise panels are used. The leading edge is 20% of 
the chord and a vorticity distribution varying linearly from 
zero to a maximum value y^ depending on the local spanwise 
load. The middle panel is 40% of the chord with constant 
vorticity y^. The trailing edge panel is 40% of the chord 
and returns the vorticity to zero linearly. For a given 
value of local loading (c*C^), the bound circulation r can 
be obtained from: 

r - 7 v c - c n 

where: V is the free stream velocity = 26.84ms”^ (as in the 

test) 

c is the local chord, meters 

is the local lift coefficient. 

T is also equal to 

\l (y 2 + Y^) < x 2 -x i> + (Y 3 +Y 2 ) ( x 3 -x 2^ + (Y 4 +Y 3 ) (x^-x-j) 

where: y^ is the value of y ^ at the ith chordwise node 

x^ is the x coordinate of the ith chordwise node. 
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Now, Y]_ = = 0, Y 2 = Y 3 = Y 2 » x 4 = 0 since the trailing 

edge is placed on the y-axis , and x-^, x 2 and x^ are respec- 
tively 1, 0.8, and 0.4 multiples of the local chord. Thus: 

^2 = TTTK ms ' 1 

Due to the linearization of c*C^, large variations were 
found in the values of the shed vorticity when computed using 
differences. A good approximation was obtained by graphi- 
cally estimating the slope of the c*C^ curve at the various 
spanwise nodes. These nodes were chosen at y = 0, 0.42, 

0.67, 0.82, 0.91, 0.96, ; 0.996, and 1.016 m, according to the 
considerations discussed in Chapter 4. 

The values of the shed vorticity y^ are obtained at 
these stations from: 

V d , „ . 

Y 1 " 1 Hy (c,C V 

The amount of shed vorticity across panel i is 



which is equal in a reduction in bound vorticity F. Thus, 
and Y 2 at each node are computed and shown in Figure 21. 
The bound vorticity is represented by the y^ values. The 
shed vorticity near the tip was initially calculated to 
satisfy Helmholtz’ law (see Figure 13), and is shown by 
point A and the dashed portion of the curve. The value was 
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32.2. Point B and the solid curve reflect a value of 45. 

The difference will be discussed later. 

Results and discussion of the validated case 

The results of the numerical modeling of the test wing 
described above are presented here. Figures 22a through 22f 
show sections of vortex sheet at various streamwise stations: 
0.35, 0.65, 0.95, 1.25, 1.55, and 2.50. The stations calcu- 
lated are 0.35, 0.65, 0.95, 1.25, 1.55, 1.88, 2.25, 2.66, 

3.1, 3.7, 4.5 and 5.5. The last three stations are in 
Region III, that is, their shape is identical to the one at 
3.1m. The plot at station 2.50 (in Figure 22f) is inter- 
polated between 2.25 and 2.66, and is used to compare with 
the available experimental results. 

First station in the sheet The trailing edge of the 
wing is at x=0, and this is where the first station should 
be. However, the linear distribution of the bound vorticity 
led to unrealistically large localized vorticity at the trail- 
ing edge, especially near the tip, and it was necessary to 
displace the first station of sheet an arbitrary 0.05 m 
downstream. The downwash angles thus calculated are more 
consistent with experience (Figure 23) . 

Tip vorticity The computed value of tip vorticity 
(point A in Figure 21) is inaccurate due to the linear 
approximations at the tip. The vorticity there is theoreti- 
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Figure 21. Piecewise linear rep: 
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Figure 22e. Section of vortex sheet at station 1.55 
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cally infinite. The results using the calculated value (A) 
are plotted using the dashed lines. The computed slope at 
the tip results in a reversal of curvature there, due to the 
reduced vorticity value and other effects discussed below. 

The tip vorticity was subsequently increased to 45 (point B) 
arbitrarily, and improved results were obtained at some of 
the stations. However, the shapes thus obtained show a 
fuller curvature and a better "spiral" shape, consistent with 
intuitive expectations of the true shape. 

Two alternatives were considered. The first one was to 
reduce the spanwise size of a few panels near the tip to im- 
prove the vorticity modeling. The result was an excessive 
proximity of the quadrature points (Imped singularities) 
and the panel nodes, resulting in numerical problems. Usage 
of double precision to reduce round off error would have 
helped only slightly, but would result in doubling the compu- 
tational resource requirements. Ideally, a closed form 
integral, albeit approximate, would probably solve this 
problem. The second was to increase this arbitrary tip 
vorticity to larger values, e.g. , 100, violating Helmholtz' 
first law, and still resulting in numerical problems similar 
to those above, now due to the large value of vorticity. 

Slope at the tip The quality of the solution was 
found to be extremely sensitive to this, variable . Originally, 
it was planned to use the chord of the tip panel to approxi- 




Figure 23. Illustration of downwash variation downstream of wing 
at the plane of symmetry 


90 


mate the tangent at the last node, and in fact, this was used 
in the preliminary testing of the method. However, flatten- 
ing and inflections in the curve near the tip resulted, as 
expected. The second attempt was to use the bisector of the 
angle between the chords of the two adjacent panels at the 
tip, and estimate the magnitude using a heuristic formula. 
While this was successful in some cases, it failed in most, 
simply because it is difficult to program logic to account 
for all possibilities of angle combinations, orientations and 
relative node positions. In addition, this method is unrelat- 
ed to the physical aspects of the problem. The method 
described in Chapter 4 was finally resorted to. The approxi- 
mation used is essentially an improvement of the tip chord 
approach. A point is chosen within the tip panel to produce, 
along with the tip node in a short chord to approximate the 
required tangent. While some numerical problems are overcome 
(namely, eliminating the need for a very small tip panel), 
others resulted. The computed induced velocity at that point 
is dependent on its position relative to the quadrature 
points. Thus, its location was varied until the least amount 
of inflections in the curves were obtained. It seems that 
the induced velocity there is still somewhat excessive. It 
also seems that the exact value of the slope (also described 
in Chapter 4) should be used for better results, that is, to 
reduce the inflections or eliminate them. 
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Comparison with test This is done in Figure 22f. 

The dashed line labelled "Vapor condensation" is at the 
center of the bands of illuminated water particles as they 
appear in the photograph in reference (29), Figure 5.4. The 
circular area labelled "Vortex core" appears as a dark circle 
in the photograph. This is the innermost part of the rolled- 
up sheet, where viscous effects are most prominent. The 
computed contour is shown with the solid line. The computa- 
tion was done starting with a flat sheet, and while itera- 
tions could have been performed to "relax" the shape to one 
where the pressure differential across the sheet vanishes, 
the closeness of the computed and experimental results sug- 
gested this was unnecessary for the purpose of this study, 
which is to test the ability of this modeling to produce 
realistic predictions of the vortex sheet. In their paper, 
Butter and Hancock (6) show the results of their relaxation 
scheme. There, the second cycle's results are close to 
those of the initial cycle, except in the two stations immedi- 
ately behind the wing, where the rolled up sheet effects are 
not included for the first, and only slightly included for 
the second. It must be noted that the only available photo- 
graph of the flow visualization is for a station far enough 
downstream, that viscous effects have had a considerable 
effect on the shaping of the sheet. Furthermore, it appears 
that, if the tip condition were modeled more accurately, a 
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better approximation would have resulted. 

Another consideration is that the test wing has a 
trapezoidal planform with considerable tip flow, not account- 
ed for in the mathematical model. 

Comparison with other computational methods Reference 

(29) shows predictions of the sheet shape at the same station, 
using a very large number of discrete vortex filaments, both 
inviscid and viscous, the results show the spiral shape is 
contained between 96 % and 101% of the half span. In other 
words, the roll-up form is too tight, compared with experi- 
ment. The calculations of these methods are based on the 
Trefftz plane scheme (two-dimensional) . The results of the 
current methods are much more realistic. 

The Wing with Deployed Flaps 

In order to assess the behavior of the present method in 
predicting vortex sheet rollup resulting from complex wing 
loading, such as for deployed partial span flaps and partially 
stalled wings, the loading of the test wing was altered by . 
increasing the lift in the inboard half, as shown in Figure 
24, by the bound vorticity curve. A slight dip in lift was 
included just outboard of the flaps to render the loading 
more realistic. The remainder is unchanged. 
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Numerical modeling 

The procedure is identical to the one just explained, 
with the exception of the selection of the number of spanwise 
nodes. They must be increased where the rates of change in 
vorticity is high, not only to model it accurately, but to 
improve the surface fit, as will become evident in the dis- 
cussion. Twelve spanwise patches were used here, as shown 
in Figure 24. 

Results and discussion 

The results are shown in Figures 25a through 25e, for 
sections of the sheet at stations 0.35, 0.65, 0.95, 1.25, and 
1.88. All but the last are very realistic, even though no 
experimental results are available for comparison. The solu- 
tion is consistent up to 1.25 half spans downstream. At 
1.88 m, the vortex sheet crosses itself. This is interpreted 
as follows: The nodal displacements are calculated properly 

for that station, however, the surface fit is done according 
to a criterion which is not necessarily consistent with the 
physics of the problem. Second derivative continuity is 
required in the fit, and this poses constraints on the poly- 
nomial degrees of freedom. Thus, the resulting solid curve 

o 

in Figure 25e is merely a C interpolant of the node points 
shown. From a physical standpoint, however, a more realistic 
constraint to impose is the slope of the tangents at the 
nodes, with values computed from the induced velocity gradient 
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Figure 24. Piecewise linear representation of typical wing 
with deployed flaps 














Figure 25e. 
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along the curve, discussed in Chapter A for the tip slope. 
If such a criterion were used at all nodes, instead of the 
tip only, the result may be similar to the dashed line in 
Figure 25e. Another way to achieve this is to increase the 
number of nodes to improve the fit, but at an increased 
computational cost. 
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CHAPTER VII: 

CONCLUSIONS AND RECOMMENDATIONS 


The objective of this study was to test the concept of 
using parametric bicubic patch surface definitions with bi- 
linear vorticity or biquadratic doublet distributions, to 
model the vortex sheet and predict its shape numerically in 
the vicinity of the wing. 

This choice stems from the need for continuity in these 
physical quantities, which is the way they occur in nature. 

To a full extent, this is true of the geometric representa- 
tion, and to a lesser extent, it is true of the singularity 
distribution. For the purpose of computing induced veloci- 
ties, it was shown that higher order of vorticity would not 
improve the truncation error. For the purpose of modeling 
the tip vorticity, the bilinear distribution, although satis- 
factory, gave rise to some difficulty. The bicubic geometric 
surface representation proved very suitable for curved sur- 
faces such as the rolled-up vortex sheet, with the surface 
fit scheme failing only for complex wings, far downstream 
beyond the region of concern. 

The author's opinion is that these difficulties can be 
surmounted by using the exact value of the velocity gradient 
at all the nodes of the sheet to obtain the tangent vectors 
in the spanwise direction. Therefore, it is suggested that 
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this effort be undertaken in a follow-up study. The. versa- 
tility of the model may also be improved by using a more 
appropriate integration scheme for handling singularities. 
Despite these approximations, the results shown in Chapter VI 
amply demonstrate the quality of the approach compared to 
some of the schemes currently used, based on discrete, piece- 
wise linear vortex filaments. For example, the Trefftz plane 
method (32) results in very tight spiral (29), and the method 
used by Hackett and Evans (14), which uses finite upstream 
filament lengths, displays excessive sensitivity to the posi- 
tioning of the filaments. The present method is sensitive to 
node positioning only near the tip and trailng edge, and the 
author's opinion is that it can be alleviated by the above 
suggested refinements. 

This method should be incorporated in a potential flow 
computational system for a variety of reasons. 

1. It is certain to improve the accuracy psie.cU.ctA.ng 
the. ootng loa.cU.ng, following one iteration, despite 
the fact that tip flow is not included. In most 
available programs, the Kutta condition is imposed 
on the assumption that the streamlines at the trail- 
ing edge coincide with some arbitrary direction, for 
example, the bisector of the trailing edge angle. 
This may be acceptable as a first estimate, giving 
downwash angles at the trailing edge to be used for 
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streamline directions for the next and probably last 
relaxation cycle. 

2. Farther from the wing, that is, near aft fuselage 
and empennage surfaces, the effect of rollup is 
necessary in many cases in attempts to pxzdlct 
stability dcaivativ Z6 , both longitudinal and 
lateral. Lateral conditions will necessitate 
modeling of both halves of the aircraft, which is 
done in potential flow solutions involving side- 
slip, yaw or roll rates. A 4 yAtcmatic, fiigonou* 
and tiational evaluation of the configuration, 
especially if unique and original, will be obtained, 
instead of attempts of intuitively adapting existing 
data which are not necessarily valid, potentially 
leading to inadequate design to be discovered only 
after expensive wind tunnel testing. 

3. kddltlo nal computational c^ofit i& ficlatlv cly 4mall, 
except in the cases requiring re-evaluation of the 
entire potential flow solution, as described in 1 
above. The solutions in 2 will not require reevalua- 
tion of the entire flow, once the wing loading has 
been determined, in general. 

4. The usage of the parametric bicubic patch in the 
manner of the present study suggests it can be used 
to represent solid boundaries as well, for the 
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purpose of solving for the singularity (source) 
strength distributions, and hence, pressure dis- 
tributions. Fewer panels will be needed than with 
flat panels, and the nodal spacing is not very 
critical due to bilinear singularities. This fea- 
ture is extremely valuable for preliminary designers, 
and its usage is being introduced (18). The com- 
pa££b£t££y of such a system with the techniques of 
the present study is evident, and incorporation of 
this wake relaxation scheme can be readily imple- 
mented. 

Additional experimental verification is desirable, 
especially for the flaps deployed case, as well as a yawed 
configuration. These seem to be natural extensions of the 
work presented here, along with the suggested refinements 
discussed earlier in this chapter and in Chapter VI. 
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